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A simplified computational method for studying the heat transfer characteristics of parallel plate ther¬ 
moacoustic heat exchangers is presented. The model integrates the thermoacoustic equations of the stan¬ 
dard linear theory into an energy balance-based numerical calculus scheme. Details of the time-averaged 
temperature and heat flux density distributions within a representative domain of the heat exchangers 
and adjoining stack are given. The effect of operation conditions and geometrical parameters on the heat 
exchanger performance is investigated and main conclusions relevant for HX design are drawn as far as 
fin length, fin spacing, blockage ratio, gas and secondary fluid-side heat transfer coefficients are con¬ 
cerned. Most relevant is that the fin length and spacing affect in conjunction the heat exchanger behavior 
and have to be simultaneously optimized to minimize thermal losses localized at the HX-stack junctions. 
Model predictions fit experimental data found in literature within 36% and 49% respectively at moderate 
and high acoustic Reynolds numbers. 
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1. Introduction 

The enhanced use of renewable energies has nowadays become 
a major concern for the world energy policy to meet the target of 
sustainable development. Thermoacoustic (TA) technology can 
play a significant role in this context because of its remarkable 
advantages over conventional energetic technologies. TA engines 
are sound-heat energy conversion devices which operate with no 
moving components, use no-polluting working gases, can be 
powered by low-grade energetic inputs (waste heat, solar energy, 
etc.) and which join construction simplicity/reliability to low 
fabrication costs. 

The physical key mechanism for these favorable technical fea¬ 
tures is the “thermoacoustic effect” where a sound wave naturally 
constrains fluid particles oscillating near solid surfaces to undergo 
properly phased pressure changes and heat transfers through a 
thermodynamic cycle [1,2]. The sound pressure-velocity phase an¬ 
gle (imposed by an acoustic network) distinguishes between trav¬ 
eling-wave and standing-wave engines, in the former the thermal 
interaction of the oscillating gas with the solid boundaries of a 
regenerator (a high heat capacity porous medium of small hydrau¬ 
lic radius) results in Stirling-type cycles with engine conversion 
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efficiencies above 40% of the Carnot efficiency [3], in the latter 
the fluid particles perform Brayton-type cycles in a stack (of larger 
pores) but "intrinsic” thermal irreversibilities limit the engine per¬ 
formances typically below 20% of Carnot’s [3]. 

For the thermoacoustic effect to be fully exploitable in power- 
production/refrigeration applications the stack/regenerator must 
work in conjunction with a couple of heat exchangers (HXs) which 
transfer heat between its edges and external heat sources and 
sinks. These components have recently been the object of many re¬ 
search studies since a significant improvement of the overall en¬ 
gine’s performance is expected to derive from their optimization. 
The finned-tube [4], and shell and tube [5] HX configurations are 
commonly used but reliable and univocal design criteria are still 
lacking. The main goal in designing efficient HXs is the achieve¬ 
ment of high transfer rates in conjunction to low acoustic dissipa¬ 
tion. Addressing this task for oscillatory fluid flows entails a series 
of new challenges compared to traditional compact HXs working in 
steady and unidirectional flows: 

(a) Longitudinal (along the x direction parallel to the particle 
oscillation) and transverse (along the y direction normal to 
the fin/plate surfaces) dimensions are bounded respectively 
by the particle displacement amplitude x^ and by the ther¬ 
mal penetration depth <5 K (= ^J2 k/<d), the distance through 
which heat can diffuse in an acoustic cycle (k and co being 
respectively the gas thermal diffusivity and the angular fre¬ 
quency of the sound wave). This imposes severe restrictions 
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Nomenclature 



a 

speed of sound (m s -1 ) 

P 

thermal expansion coefficient (K -1 ) 

A 

surface area (m 2 ) 

y 

ratio of isobaric to isochoric specific heat 

Cp 

isobaric specific heat (J kg -1 1C 1 ) 

s 

penetration depth, m 

c F 

specific heat of the secondary fluid (J kg -1 K ’) 

n 

shear viscosity coefficient (N m -2 s -1 ) 

DR 

drive ratio (=Pa/Po) 

0 

TWo (K) 

e 

energy flux density (W m -2 ) 

AT HX s 

stack-HX T difference (K) 

f 

spatially averaged h function 

AT H xg 

gas-HX T difference (K) 

h 

thermoviscous function, specific enthalpy (J m 3 ), heat 

K 

gas thermal diffusivity (m 2 s -1 ) 


transfer coefficient (W m -2 K 1 ) 

X 

wavelength of the sound wave (m) 

j 

imaginary unit 

p 

h F P T /mc F 

jc 

Colburn-j factor 

V 

kinematic viscosity (m 2 s -1 ) 

k 

wave number (m -1 ), thermoviscous function (m) 


generic acoustic variable 

I< 

thermal conductivity (W m -1 K ') 

n 

thermoviscous function 

l 

half plate/fin thickness (m) 

p 

density of the gas (kg m -3 ) 

L 

heat exchanger or stack length (m) 

<z> 

viscosity stress tensor (N m -2 s -2 ) 

rh 

secondary fluid mass flux (kgs -1 ) 

X 

fit constant of experimental values of I( s , K c and I< H (K -1 ) 

Mol 

gas molecular weight (kg kmol ') 

CO 

angular frequency (rad s -1 ) 

Nu 

Nusselt number 

C2 

blockage ratio 

P 

acoustic pressure (Pa) 



P 

pressure (Pa) 

Subscripts 

Pr 

Prandtl number 

0 

mean, time averaged 

P T 

internal perimeter of a tube section (m) 

1 

first order acoustic variable 

<? 

time-averaged heat flux density (W m 2 ) 

A 

acoustic amplitude at a pressure antinode 

Q 

heat load (W) 

C 

referred to the cold HX 

Q* 

heat flux at the fin surface (W) 

F 

secondary fluid 

R 

gas universal constant (J K 1 kmol -1 ) 

H 

referred to the hot HX 

R' 

R Mol -1 (kj K -1 kg -1 ) 

in 

evaluated at the inlet of the HX 

t 

time (s) 

p 

isobaric 

T 

temperature (K) 

res 

resonator 

u 

specific internal energy (J m -3 ) 

ref 

reference value 

u 

fin-to-reservoir unitary global conductance 

s 

stack 


(W m -2 K -1 ) 

solid 

solid 

V 

acoustic particle velocity (m s -1 ) 

X 

longitudinal, x-component 

X 

axial coordinate/direction (m) 

y 

transversal, y-component 

X-S 

stack center coordinate (m) 

Z 

along the z direction 

y 

transverse coordinate/direction (m) 

K 

thermal 

yo 

half plate/fin spacing (m) 

V 

viscous 

Greek symbols 

Superscripts 

a 

fit constant of experimental values of I< s , I<c and I( H 

in 

evaluated at the inlet of the HX 


(W m -1 K -1 ) 

* 

evaluated at the fin surface 


on the heat transfer area and on the flow resistance 
minimization. 

(b) The junctions between the stack/regenerator and the HXs 
are interested by complex temperature and flow patterns 
due to cross-sectional abrupt changes and entrance effects 
[6,7], These phenomena may strongly affect the actual heat 
transfer and power dissipation rates. 

(c) Unambiguous and reliable heat transfer correlations laws for 
the oscillatory flow regime in the sub-kHz domain are not 
yet available. Conventional steady flow correlations are 
commonly applied for the estimation of the heat transfer 
coefficients previous "ad-hoc” modifications [8,9]. 

(d) Large temperature differences are imposed by the hot and 
cold HXs to the ends of the stack/regenerator (typically long 
only several centimeters). These high temperature gradients 
immediately results in performance degradation (especially 
in TA refrigerators) due to the associated increased thermal 
irreversibilities. Temperature differences between the HXs 
and the stack/regenerator should be minimized preserving 
however an efficient inter-element heat exchange. 


These issues cannot be faced by the standard linear theory [1,2] 
owing its ID character and being based on the so called “mean 
field approximation” [10]. Thus 2D or 3D numerical and analytical 
models able to account for the transverse heat transfer phenomena 
taking place in the HXs are currently developed [11-15]. In this pa¬ 
per the results of a 2D numerical model for the study of the ther¬ 
mal performance of parallel plate thermoacoustic HXs coupled to 
a stack in a simplified configuration are presented. The model de¬ 
scribes resonator acoustics through an idealized ID sound wave 
while performs detailed 2D computations of the time-averaged 
temperature and energy flux distributions within a representative 
domain of the stack and HXs. These last are generated by an energy 
balance-based numerical calculus scheme which implements the 
energy transport variables of the classical linear thermoacoustic 
theory. The model constitutes an extended version of the one pre¬ 
sented in [8,16] (simulating a thermally isolated stack) by inclu¬ 
sion of hydrodynamic heat transport along the transverse y 
direction, temperature dependent thermophysical gas/solid 
parameters, viscous dissipation and adjacent HXs interacting with 
“hot" and “cold” thermal reservoirs. The model is validated by 
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comparing numerical predictions against experimental data found 
in literature. The effect of operating conditions and geometrical 
parameters on the heat transfer characteristics are investigated. 
Information on the optimal settings which maximize the HX per¬ 
formance is derived. 

2. Formulation 

The thermoacoustic model system considered in this work is a 
stack of parallel plates, of length L s , sandwiched between “hot” 
and “cold” parallel-fin HXs respectively of length l H and L c . This 
assembly is located at a distance x s from the center of a half-wave¬ 
length (A/2) gas filled resonator sustaining a standing acoustic 
wave as shown in Fig. 1. Stack plates and HX fins are modeled with 
the same spacing (2y 0 ) and thickness (21) so the three structures 
result characterized by the same blockage ratio Q = A s IA res = 1 / 
(1 + 1/yo), where A res is the cross sectional area of the resonator 
and A s is the cross sectional area of the stack/HXs open to gas flow. 
This choice, while allowing for a simple geometry, should lead to a 
reduction of the detriment effects discussed in Section 1 at point b, 
being the flow and energy path almost not interrupted at the junc¬ 
tions. The implicit assumption is that an efficient HX is being mod¬ 
eled. The HX fins are set close to the ends of the stack plates but are 
not in touch with them in order to reduce adverse heat conduction 
leaks from hot to cold HX. To include in the calculation the effect of 
the heat transfer coefficient from the secondary fluid-side, the HXs 
are modeled to thermally interact with two thermal reservoirs, the 
“hot" reservoir at temperature T H and the “cold” reservoir at tem¬ 
perature T c . Since we are considering a thermoacoustic device 
working in the refrigeration mode the former acts as a heat sink 
while the latter acts as a heat source. 

In the rest of this section the basic energy balance equations on 
which the model relies upon are defined together with the expres¬ 
sions of the energy transport variables entering in them. Prelimi¬ 
narily, the main assumptions involved in the model are discussed. 

2.J. Basic approximations and simplifying conditions 

The formulation relies on the following assumptions: 

• The problem is two-dimensional; 

• The acoustic approximation is valid (low Mach number regime) 
and any acoustic variable c can expresses in complex notation 
by the conventional first-order expansion {(x.y.t) = £ 0 (x,y)+ 
Reji^x.y)^} where t is the time, j the imaginary unit, Re{} 
signifies the real part and where the mean value £ 0 is real but 
the amplitude is complex to account for both magnitude 
and phase of the oscillation. 

• The working fluid is a Newtonian gas with viscosity obeying the 
ideal gas state equation P 0 /po = R'T 0 where P 0 is the static pres¬ 
sure, p 0 the static density, T 0 the mean (time-averaged) temper¬ 



Fig. 1 . Schematic illustration of the modeled thermoacoustic refrigerator. 


ature and R' = R/Mol (R being the universal constant of gases and 
Mol the molecular weight of the gas under study). Temperature 
dependence of the gas thermophysical properties p (thermal 
expansion coefficient), 17 (shear viscosity coefficient) and K 
(thermal conductivity coefficient) are accounted for by the 
relations 


"-h "-"-’ijtf- K=K «&T 

J 7 re/and K re f being the values of i; and K at the reference temper¬ 
ature T ref . 

• The stack is significantly shorter than the acoustic wavelength 
and is not intrusive. This allows the sound field to be retained 
as spatially uniform over the HXs/stack assembly and to be 
approximated in the neighbor (but also inside the stack as far 
as the pressure is concerned) by a ID loseless standing wave: 

Pi = Pa sin (kx s ) = p 0 , v x , <=j^~ cos(kx s ) = jv, (2) 

Po a 


p, being the amplitude of the dynamic pressure, P A the ampli¬ 
tude of the dynamic pressure at a pressure antinode, v x \ the 
amplitude of the longitudinal particle velocity, k the wave num¬ 
ber (k = 27t/A) and a the sound velocity. 

• The standard thermoacoustic equations are considered for the 
amplitudes of the oscillating temperature Ti, longitudinal veloc¬ 
ity zy,, transverse velocity gradient 6 v xl / 8 y transverse velocity 
v y \ and transverse velocity gradient 8 v yl /Sy, inside a gas pore. 
If the specific heats of the plate/fin materials are notably greater 
than the isobaric specific heat of the gas c p (negligible temper¬ 
ature oscillations in the solid), the first three equations read as 
[ 1 . 2 ]: 


t 1 , , 1 dp , <9To 

1 - p 0 c P K ^ 1 p 0 co 2 (1-Pr) dx dx K 

- Pr( 1 - h v )] 


(3) 


Vxl 


j dp, M h , dv A _ 2 dpi h 

cop, dx nv) ’ dy cop 0 8 2 v dx kv 


where 


= cosh[(l +j)y/S K ] 

K cosh[(l +j)y 0 /S K ] ’ 


= cosh[(l +j)y/8 v ] 
’’ cosh[(l +j)y 0 /$v] 


and 


(4) 


5 k sinh[(l +j)y/S K ] 5 V sinh[(l +j)y/8 v ] 

K 1 +j cosh[(l +j)y 0 /d K ] ’ v 1 +j cosh[(l +j)y 0 /5 v ] 

and where Pr is the Prandtl number, <5 V = sjlv/co the viscous pen¬ 
etration depth (v being the kinematic viscosity) and where y = 0 is 
fixed at the center of the fluid gap. 


The analytical expression of zy, can be derived from the linearized 
equation of continuity 


J " Pl + It Vxy ) + dy (Po Vn ) = ° 


(5) 


Neglecting the term (0po/8y)zyi in comparison to the term (Spo/ 
6 x)zy,, (Vyi being expected to be of order c5 v /A smaller than v x f] 
and taking into account the relations 


1 9p 0 
Po dx 


Tg dx 


Pi =-PoPTi+-^P 


(y being the ratio of isobaric to isochoric specific heat) Eq. (5) can be 
put into the form 
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dVy\ 

dy 


=- 


dv x \ 

dx 



( 6 ) 


dv n . co f [1 + ( y - 1M1 - h r ) - [1 + (y - l)h K ](l -/ T ) 

(1 ~fv) 

7,.(1 - h K ) -f K ( 1 - h,) + ( h K - h v ) 


u yc _j cu \[ 

dy J p 0 a 2 \ 


P i 


+] 


p 0 co{\-Pr) _ 


(1 -fv) 


9To _ 
dx dx [ 


whose integration between 0 and y provides 


„ , ® /[i + (7-i)/ K ](y-fcv)-[y+(y-i)^](^-/v) 

y yi =J ^2 


p 0 a 
+J 


Pi 


p 0 co(l-Pr) 


(1 -/v) 

7v(y - k K ) - f K (y - fey) + (fc K - fey) ] dT 0 dp, 
(1 -fv) j dx dx 


( 8 ) 


where 


_ tanh[(l + j)y 0 /d K ] _ tanh[(l +j)y 0 /d v \ 
jK [(i +i)y 0 /s K ] ■ Jv [(1 +j)y 0 /Sy\ 

In Fig. 2 the profile of the real and imaginary parts of Eq. (8) to¬ 
gether with its module is reported as a function of the transverse 
coordinate y for parameters corresponding to run 1 of Table 1. 
The transverse velocity vanishes both at mid channel y = 0 (for sym¬ 
metry) and at the wall surface y = y 0 . 

• The longitudinal pressure gradient dp,/dx along the stack/HXs 
assembly is calculated by imposing continuity of volume flow 
rate at the entrance of the HXs 


Eqs. (3) and (4) can now be substituted into Eq. (6) for Ti and v x \ 
while the axial velocity gradient can be deduced from the wave 
equation (Eq. (A19) of Ref. [1]). Making substitutions and rearrang¬ 
ing (by also using the thermodynamic relation c p = /3a 2 /(y - 1)) we 
obtain 


dPi _ Po® 
dx (l-7 v )fl 0 


( 9 ) 


We remark as the validity of the above derived equations could be 
invalidated at the pore ends of the stack or at the HX fins where the 


Table 1 

Parameters of selected simulations. In all runs test-gas =He, P 0 = 101325 Pa, mean temperature = 300 K, resonance frequency = 200 Hz, = 100 W m 2 K 1 , A= 10.08 m, 
S K = 5.352 x 10- 4 m, T H = 300 K, T c = 297 K, L s = 0.07 m. 


Run 

yo 18k 

1I&K 

XsP 

Pa/Po (%) 

Lc/Ts 

LhILs 

Uc/Uref 

UnlUref 

1 

3 

0.93 

0.11 

4.93 

0.17 

0.18 

1 

1 

2 

1.5 

0.47 

0.11 

4.93 

0.11 

0.11 

0.1 

0.1 

3 

1.5 

0.47 

0.11 

0.49-5.43 

0.11 

0.11 

0.1 

0.1 

4 

1.5 

0.47 

0.11 

0.49-5.43 

0.11 

0.11 

1 

1 

5 

1.5 

0.47 

0.11 

0.49-5.43 

0.11 

0.11 

5 

5 

6 

1.5 

0.47 

0.11 

3.95 

0.11 

0.11 

0.1-20 

0.1-20 

7 

1.5 

0.47 

0.11 

4.44 

0.11 

0.11 

0.1-20 

0.1-20 

8 

1.5 

0.47 

0.11 

4.93 

0.11 

0.11 

0.1-20 

0.1-20 

9 

1.5 

0.47 

0.11 

5.43 

0.11 

0.11 

0.1-20 

0.1-20 

10 

1.5 

0.47 

0.11 

2.47 

0.11 

0.11 

0.1-20 

0.1-20 

11 

1.5 

0.47 

0.11 

3.45 

0.11 

0.11 

0.1-20 

0.1-20 

12 

3 

0.93 

0.11 

0.49-5.92 

0.11 

0.11 

10 

10 

13 

1.5 

0.47 

0.11 

0.49-5.92 

0.11 

0.11 

20 

20 

14 

1.5 

0.47 

0.11 

0.49-7.90 

0.11 

0.11 

20 

20 

15 

1.5 

0.47 

0.11 

2.96 

0.0143-0.343 

0.18 

30 

30 

16 

1.5 

0.47 

0.11 

4.93 

0.0143-0.343 

0.18 

30 

30 

17 

1.5 

0.47 

0.11 

6.91 

0.0143-0.343 

0.18 

30 

30 

18 

0.5-3 

0.16-0.93 

0.11 

4.93 

0.0143 

0.18 

30 

30 

19 

0.5-3 

0.16-0.93 

0.11 

4.93 

0.036 

0.18 

30 

30 

20 

0.5-3 

0.16-0.93 

0.11 

4.93 

0.071 

0.18 

30 

30 

21 

0.5-3 

0.16-0.93 

0.11 

4.93 

0.126 

0.18 

30 

30 

22 

0.5-3 

0.16-0.93 

0.11 

4.93 

0.180 

0.18 

30 

30 

23 

0.5-3 

0.16-0.93 

0.11 

4.93 

0.343 

0.18 

30 

30 

24 

2 

0.8 

0.11 

1.63 

0.14 

0.14 

20 

20 

25 

2 

0.8 

0.11 

6.91 

0.14 

0.14 

20 

20 

26 

0.5-3 

0.75 

0.11 

4.93 

0.036 

0.18 

30 

30 

27 

0.5-3 

0.75 

0.11 

4.93 

0.018 

0.18 

30 

30 

28 

0.5-3 

0.75 

0.11 

4.93 

0.343 

0.18 

30 

30 

29 

1.5 

0.47 

0.09 

7.9-8.9 

0.29 

0.29 

20 

20 

30 

1.5 

0.47 

0.09-0.08 

8.9 

0.29 

0.29 

20 

20 

31 

1.5 

0.47 

0.07 

9.9 

0.29 

0.29 

20 

20 

32 

1.5 

0.47 

0.11 

0.99-7.89 

0.11 

0.11 

20 

20 

33 

1.5 

0.47 

0.1-0.09 

7.89 

0.11 

0.11 

20 

20 

34 

2.0 

0.8 

0.11 

0.69-8.0 

0.14 

0.14 

20 

20 

35 

2.0 

0.8 

0 . 11 - 0.1 

8.0 

0.14 

0.14 

20 

20 

36 

2.0 

0.8 

0.09-0.07 

8.0 

0.43 

0.43 

20 

20 

37 

2.0 

0.8 

0.07 

8.4-10.8 

0.43 

0.43 

20 

20 
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gas mean temperature is expected to become strongly y-dependent 
to allow net heat exchange with the solid surface. Eq. (3), instead, 
was derived as solution of the heat transport equation involving a 
y-independent T 0 . Analogously, invoking validity of Eqs. (4), (7) 
and (8) at the fin outer ends (facing the duct) corresponds to neglect 
any entrance effect which may be the source of minor losses. Re¬ 
sults of the present study have thus to be appreciated as a function 
of these major limitations. 

2.2. Energy balance equation and energy fluxes in the gas 

Under the simplifying conditions outlined in the previous sub¬ 
section the gas flow dynamics remains completely described 
through Eqs. 2, 4, and 9. Meanwhile, the temperature field is gov¬ 
erned by the general equation of energy conservation [17] 

t L(lpv 2 +pu^=-Ve (10) 


where tilde indicates complex conjugation. Substituting now Eqs. 
(2) and (9) into Eqs. (3), (4), (7) and (8) and these last in Eq. (14), 
the following expression is found, at second order in the acoustic 
oscillation amplitude, for e x : 


e * = 2fl lm 


n,n 2 

(1 -fv) 


PoVo- 


Po°p 


dT 0 


+ 


11 CO 


3 Qp 0 a 2 

OTo 


Re 


(1 -fv) 


2co£2 2 (l-Pr)|l-/ v | 2 dx 
ijP 

Po v 0 + —5- - --T 

312 2 (1-Pr)|l-/ V | 2 


Im 


rjco 




Im 


>1P 


dT o 


<5 2 fl 2 (l-Pr)|l-/ v | 2 dx 


k v n 2 

(1 -fv) 


lm{k r n 4 }vl-K 


PoVo 

dTo 

dx 


|n,n 2 }^ 


(15) 


where 


17, = (1 - h K ) 


where u is the specific internal energy and where e, the energy flux 
density, is defined as 

'\ 


e = p\l - v + h) - v 4> - KVT 


01) 


n 2 = (i — h v ) 

n3 = [±±fti i)/*](y - fc v) - [y + (r - i)k>c](i -fv) 


h being the specific enthalpy and <P the viscous stress tensor. Time 
averaging Eq. (10) over an acoustic cycle leads to 

V £ = 0 (12) 

where over-bar means time-averaged. Eq. (12) is applied in this 
work to each computational cell of the simulation domain to im¬ 
pose local energy balance in the gas. To accomplish this, a finite dif¬ 
ference technique is employed where the quantitative results of 
standard linear theory are considered for x and y the components 
of the time-averaged energy flux density. These last are derived in 
the following subsections. 


2.2.1. Energy flux density along the x direction 

The x component of the time-averaged enthalpy flux density is 
[17] 


V co f 

ex = 2^V 


c-Iti/w /-j \ c-2n/(o 

jPV* (2 V 2 + h J dt - J ( v x o xx + VyT xy )dt 


-I< 


L 


lulus rsj 

d~x dt 


(13) 


where the components of the viscous stress tensor are 


, dv x dv. 


CT XX =-^ 2—^-—f 


dx dy 


2 dv v 


3 ' dy 


T x y 


= Ty X = 11 




>1 


dv x 

dy 


Oyy 




4 

3 fl 


dVy 


the x velocity gradients having being neglected because of order <5 V / 
k smaller than they velocity gradients. Although the second integral 
on the right hand side of Eq. (13) is generally neglected because 
including terms of second order in velocity, it is retained in this 
work because we are trying to identify the impact of loss mecha¬ 
nisms on the performance of the two HXs. Making use of the com¬ 
plex notation and retaining only terms up to second order Eq. (13) 
can be written as: 


1 


1 


dVy 


t Po c pR e {7i Pu} + q'/Re< —df— v x \ \ --iiRe\——Vyi \ -I<— 


dy 


1 


dv x 


dy 


r dT o 


dx 

(14) 


n 4 = 

tt 5 = 

7T 6 = 


ffy - k K ) -f K (y - k v ) + (k K - k v ) 

(1 ~fv) 

[i + (y - iMi - h») - [i + (7 -1)^(1 -fv) 

(1 -fv) 

f>( 1 - K) -AO - K) + (h K - h v ) 

(1 -fv) 


2.2.2. Energy flux density along the y direction 

The time averaged hydrodynamic enthalpy flux along the trans¬ 
verse direction is 
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Proceeding analogously as done before for e x we obtain 
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which, by substitution of the explicit expression of the acoustic 
variables, becomes 
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Note as in the expressions derived for e x and e y the dependence on 
the stack position inside the resonator is provided by p 0 and v 0 . In 
the rest of this study these energy flux densities are identified to 
heat flux densities ( e x = q x , e y = q y ) because in the short stack 
approximation the work flux contribution to the total energy flux 
is small [1] (the over-bars having been neglected to simplify the 
notation). 


2.3. Energy balance equation and energy fluxes in the solid 


The energy equations in the gas are coupled to the ones in the 
solid through heat exchange at the solid surfaces. The analogue 
of Eq. (12) for the solid portions is 

V • q = 0 (18) 

where 


n K 9T ° 

Vx -~ Ks W 

and 


dT 0 


dT 0 


fl* = ~Kh —, flx = -l<c 


dx 


dx 


(19) 


fly = -I(< 


dTo 
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fly = — I<H 


dTo 
dy ’ 


fly = -K t 


dTo 

dy 


( 20 ) 


of the tube, h F the convective coefficient fluid-wall, rh the mass flux 
and c F the fluid specific heat). The total heat flux exchanged by the 
secondary fluid is 

Q. = [ h F P T 0(z)dz = mc F [ 1 - exp {-/.iL z )]6 in 
Jo 

So, the average heat flux density normal to the internal surface of 
the flat tube can be expressed in the form 

fly = = U(Tc — To) ( 21 ) 

where 

U = gj[l-exp ( -^ z )] 

plays the role of a fin-to-reservoir unitary global conductance if T in 
is identified with T c or T H (in the hypothesis that the secondary fluid 
make perfect thermal contact with the thermal reservoirs). Eq. (21) 
is used in this work to couple the heat flux thermoacoustically gen¬ 
erated in the stack to the heat carried out by the secondary fluid. 

2.5. Numerical scheme 


K s being is the thermal conductivity of the plate material and I< H and 
K c the thermal conductivities respectively of the "hot” and “cold” 
HXs. The temperature dependence of these parameters is described 
by empirical laws of the form I< so i id = K re f + a exp(±xTo)> with K re f, a 
and x determined from the fit of experimental data found in 
literature. 

2.4. Energy balance in the heat exchangers 

The cooling power Q c produced thermoacoustically in the 
stack is supposed to be extracted through the cold HX from a cold 
thermal reservoir at temperature T c . Analogously, the sum of Q c 
and of the heat produced by viscous dissipation over the surfaces 
of the stack/HXs assembly is rejected through the hot HX to a hot 
thermal reservoir at temperature T H . The modeled HXs are of the 
flat tube type (as those considered in [18]) with secondary fluids 
flowing inside the tubes along the z direction (see Fig. 3). This 
configuration leads to a more simplified modeling compared to 
other geometries but introduces in the problem a z dependence 
not included in our 2-D (x, y) treatment. To overcome this limita¬ 
tion we assume for the inner surface of the HXs tubing a constant 
temperature wall boundary condition (compatible with the high 
thermal conductivities of the HX material) but substitute for the 
local heat flux density normal to the inner tubing surface q y (z) a 
proper mean value obtained by averaging q y (z) over the length 
L z of the tube. We start from the well-known expression of the 
temperature variation T^z) of a fluid along a channel with con¬ 
stant wall temperature T 0 

d(Z) = Bin exp {-HZ) 

where 0(z) = T F (z) - T 0 , 9 in = V F - T 0 and ft = h F P T /mc F (T F being 
the inlet fluid temperature, P T the internal cross section perimeter 




The periodicity of the stack and HXs structures along y allows 
calculations to be performed in a single channel of the stack/HXs 
assembly and in the couple of plates/fins enclosing it. The compu¬ 
tational domain is further reduced by symmetry from half a gas 
duct to half a plate/fin as indicated in Fig. 4 by the light gray area 
together with the coordinate system used. The axis parallel to the 
plates is the x axis; x = 0 is chosen to be the beginning of the cold 
HX on the left. They axis is perpendicular to the stack-plates; y = 0 
is chosen to be the midpoint between the two adjacent plates or 
fins. The calculation of the steady-state two-dimensional time- 
averaged temperature distribution is performed using a finite 
difference methodology where temperature spatial gradients are 
discretized using first order nodal temperature differences. To this 
end, the computational domain is subdivided using a rectangular 
grid. In the x direction the computation mesh size is typically 
0.005 I s while in the y direction the computation mesh size is typ¬ 
ically 0.02 y 0 . The following boundary conditions are imposed: 


• Symmetry on nodal lines at y = 0 and y = y 0 + / 



9To 

dy 


y=ya+i 


= 0 


( 22 ) 


• continuity of temperature and transverse heat fluxes at the gas- 
solid interfaces (y = yo) 



Fig. 3. Schematics of the HXs and stack configuration. 


Fig. 4. Magnified view of the computational domain (gray area). 
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(23) 

where subscript “solid” refers in general both to the stack plate 
and the HX fins. 

• continuity of transverse heat fluxes at the inner surface of the 
HX tubing: 


per (I( c = I< H = 401 W itT 1 K _1 at T re /= 300 I<). In all cases the tem¬ 
perature of the cold and hot reservoirs are fixed respectively to 
T c = 297 K and T H = 300 K. The dynamic Reynolds number (t» xl t> v / 
v) and Mach number (P A lp 0 a 2 ) resulted always smaller respectively 
than 500 and 0.1, the values over which turbulence an non-linear 
behavior are expected to affect fluid dynamics [1,21]. 

The cooling power, Q _ c , is calculated by integrating q y over the 
inner tubing surface of the cold HX 


*c(^ W = Uc(Tc ~ To) - I<H(^) tube = U h (T h - T 0 ) 

(24) 

• Vanishing heat flux at the solid HX fin terminations 
Cyo ^ y ^ To + 0 facing the duct 



dTo 

dx 


= 0 


X He I l. S I 1-H ) 


(25) 


• Vanishing flux at the HX pore ends (0 ^ y <yo) facing the duct 
(gas oscillations outside the HXs are adiabatic and thus the 
thermoacoustic effect vanishes) 

q x (x = 0) = q x (x = L c + L s + L h ) = 0 (26) 

As for the separation gap between the HX fin and the stack edge, it 
is set equal to S K and, as a further simplifying hypothesis, the slice 
of gas enclosed in it is assumed at rest. 


Qc = f k U c [Tc - T 0 (x)]dx (27) 

Jo 


with L z arbitrarily set to 1 m. This quantity differs from the heat flux 
entering/leaving the outer fin surface 



by the heat flux transferred through the vertical fin surface facing 
the stack edge. Analogous quantities are calculated at the hot HX 
fin. At the operating conditions considered in this work Q c resulted 
smaller than Qf typically by 3%, evidencing as the stack perfor¬ 
mance is degraded by the heat conduction losses at the stack-HX 
junctions. 

3.1. Effect of temperature difference 


Substituting Eqs. (15) and (17) in Eqs. 12, 19, and 20 in Eq. (18) 
to impose local energy balance in each cell of the computational 
grid, a system of quadratic algebraic equations in the unknown 
variable T 0 is generated. The system is solved by a code developed 
by the author in FORTRAN-90 language which executes the recur¬ 
sive Newton-Raphson method [19]. At each iteration the system of 
linear algebraic equations providing the temperature corrections 
for the next step is solved using a LU decomposition with partial 
pivoting and row interchanges matrix factorization routine. The 
latter is taken from the LAPACK library routines available online 
at [20], where details about accuracy, computation cost, etc. can 
be found. Once the time-averaged temperature distribution is 
known, it can be substituted in Eqs. (19) and (20) to determine 
the energy flux distributions along the x and y directions both in 
the gas and in the plate/fins. In addition, it can be used to calculate 
other variables of interest such as the energy dissipation, the en¬ 
tropy generation, etc. 

Compared to other approaches currently used to simulate 
thermoacoustic devices the proposed model has the advantage of 
a reduced computational cost since it directly deals with time- 
averaged quantities and does not require integration of the conser¬ 
vation laws of mass, momentum and energy over the large number 
of acoustic cycles driving the system towards steady conditions. In 
Section 3.5 the code is validated by comparing numerical predic¬ 
tions against experimental data found in literature. 

3. Results and discussion 

The performance of the coupled stack-HXs system is investi¬ 
gated at varying both operating conditions (DR = P A IP 0 , x s , U H , U c ) 
and stack/fin geometry (y 0 , /, L H , L c ). The parameters of different 
numerical simulations (runs) are listed in Table 1. in all tests he¬ 
lium at atmospheric pressure and at a mean temperature of 
300 K is assumed as working fluid (v = 122 x 1CT 6 m 2 s _1 , 
K= 0.152 W nr 1 r 1 , a = 1008 m s“\ Pr= 0.68 at T ref = 300 K). It is 
considered enclosed in a half-wavelength resonator (acoustic 
wavelength X = 5.04 m) having a resonance frequency of 200 Hz; 
The stack is 0.07 m long and is constituted by stainless steel 
(K s = 14.9 W m _1 K _1 at T re f= 300 K) while the HXs material is cop¬ 


Fig. 5 illustrates the longitudinal temperature profile evaluated 
in correspondence of the centerline (y = yo + 1) of the plate and of 
both fins for a drive ratio DR (=P A /P 0 ) = 4.93% (run 2). The temper¬ 
ature distribution is linear in the stack (except at the stack edges) 
becoming flat in the fins, consistently with their high thermal con¬ 
ductivity (near isothermal fins). Temperature discontinuities A Tfjxs 
are observed between the fins and the stack edges (with the stack 
temperature higher than the cold fin temperature) in accordance 
with experiments [22]. At low U values (U < 100) AT HX s increases 
at increasing drive ratios as shown in the insert of the same figure 
(run 3). This behavior can be interpreted observing that in the pres¬ 
ent model the HXs are thermally interacting with hot and cold res¬ 
ervoirs. So, depending on the magnitude of the fin-to-reservoir 
overall heat transfer coefficient U, the cold fin must cool down be¬ 
low T c and the hot fin must heat up above T H until the heat flux ex¬ 
tracted from the cold reservoir and the one delivered to the hot 
reservoir balance the heat flux thermoacousticaily pumped along 



Fig. 5. Centerline plate and fins temperature as a function of coordinate x. In the 
insert the temperature discontinuities between the cold fin and the stack edge is 
reported as a function of drive ratio at selected U values. 
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the stack. At low U values this causes an increase of the mean lon¬ 
gitudinal temperature gradient across the HXs-stack assembly and 
thus of the temperature discontinuities observed at the HXs-stack 
junctions. The same figure shows as the effect becomes less 
marked at increasing U values (run 4) till to invert its trend (A T HX s 
decreases with DR) for U> 100 (run 5). In this latter case, in fact, 
the temperatures of the hot and cold HXs are “anchored” respec¬ 
tively to T h and T c and the enhanced hydrodynamic enthalpy flux 
activated by increasing DR are compatible with lower axial tem¬ 
perature gradients. This is clearly shown in Fig. 6 where the cooling 
load results to be a decreasing (linear) function of the ratio be¬ 
tween the temperature of the cold HX and the temperature of 
the hot HX (runs 6-9). Fig. 5 also reveals that at high U values 
AThxs vanishes in correspondence of a given drive ratio and then 
becomes negative (stack temperature lower than fin temperature). 
This circumstance, also experimentally evidenced by Mozurkewich 
[9], could allow for minimization of the detriment heat leakage 
phenomenon from hot to cold HX discussed in the previous sub¬ 
section. The U magnitude, therefore, directly affects the tempera¬ 
ture difference across the stack and this in turn influences the axial 
hydrodynamic enthalpy flux, i.e. ultimately, the cooling load. High 
U values are desirable to improve the performance of the device, as 
evidenced in Fig. 7 where the growth of the cooling load with U at 
selected DR is observable (runs 7, 9, 10, 11). 

The existence of inter-element heat transfer under zero temper¬ 
ature difference is not surprising because the driving temperature 
difference for heat extraction/immission at the HXs is the one be¬ 
tween the HXs and the adjacent gas. High values of this transverse 
temperature gradient should be avoided, however, because they 
could enhance the thermal irreversibility associated to the heat 
transfer process and, in addition, could determine a steep axial 
temperature gradient in the gas which, as argued before, reduces 
the useful cooling load. In Fig. 8 the mean gas-fin temperature dif¬ 
ference A T HXg = ( T HX c) - ( T gC ) at the cold HX is plotted as a function 
of DR at selected y 0 values (runs 12,13). Quantities (T HX c) and (Tgc) 
represent the spatially averaged temperatures respectively of the 
cold HX (evaluated at the gas-solid interface y=y 0 ) and of the 
gas (evaluated at the centerline of the gas channel y = 0): 

(Thxc) = 1 I'' T 0 (x,y 0 )dx, (T gC > = i f T„(x, 0 )dx (29) 

TC Jo L C Jo 

The plot clearly shows the existence of two regimes: one at low DR 
where A T HXg exhibits a squared dependence on DR (AT HXg oc DR 2 ) 



Fig. 6. The cooling load as a function of the ratio between the temperature of the 
cold HX and the temperature of the hot HX at selected drive ratios. Straight lines 
result from a linear fit of the simulation points. 



Fig. 7. The cooling load as a function of the fin-to-reservoir unitary global 
conductance at selected drive ratios. Continuous lines are guides for the eye. 


and one at high DR where A T HXg depends linearly on DR 
{AT HXg oc DR). This result seems to conflict with the statements of 
Swift [23] which predicts rather a linear behavior at low DR and a 
squared-one at high DR. The disagreement can be explained analyz¬ 
ing the Swift’s predictive formula based of the "thermal boundary- 
layer” theory [23]: 


AT H xg = 


Q-C^k 

~KA 


(30) 


A ( =L C L 2 ) being the heat transfer surface area. At low DR (for 
2x, -c L c ) it is expected Q c oc DR 2 so if 2x^ (cc DR) is substituted for 
L c in A (as done in [23]) the linear behavior is found while if A is 
set always constant to L C L Z (as done in the present work) the qua¬ 
dratic behavior is obtained. At high DR (for 2x^ » L c ) if the func¬ 
tional relationship Q c oc DR 2 is retained (as done in [23]) and A is 
set equal to L C L Z the quadratic behavior is found but if the depen¬ 
dence of Qc on DR is not quadratic a deviation from the attended 
behavior should be expected. This is indeed the case of the present 
study where simulations predict in the high drive ratio regime 



Fig. 8. The mean temperature difference between the gas and cold HX fin as a 
function of DR for two operating conditions. Dotted and continuous lines result 
respectively from a quadratic (low DR) and linear (high DR) fit of the simulation 
points. 
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(when 2x\ » L c ) a roughly linear dependence of Q c on DR (see next 
sub-section). 

An additional important information that can be drawn from 
Fig. 8 is that the choice y 0 l <5^ = 1.5 in conjunction with 
U = 2000 W ur 2 K _1 entails lower transverse temperature gradi¬ 
ents and higher cooling loads (results not shown) compared to 
the other. So for a given HX length, the benefit of reduced thermal 
irreversibilities could be met in conjunction with enhanced cooling 
loads for suitable HX configurations. 

3.2. Effect of drive ratio 

To investigate the effect of drive ratio on the cooling load Q c 
numerical simulations are performed at fixed fin lengths L c and 
L h . Results are shown in Fig. 9 where the mean heat loads per unit 
fin area q c and q H (and their difference in the insert) are plotted as 
a function of DR (run 15). The simulations reveal that at low drive 
ratios (2x, -c L c , L H ) both q c and q H increase proportionally to DR 2 
as stated by the linear theory. At high drive ratios (2xi » L c , L H ), 
however, the behavior of the heat loads differentiates: while q H 
continues to grow faster than linearly (even if with a reduced 
strength), q c exhibits a roughly linear dependence on DR. As al¬ 
ready pointed out by Swift [23], a candidate cause for the appear¬ 
ance of the linear behavior could be the inefficiency of the HXs in 
transferring the full thermal load at high pressure amplitudes. This 
should happen when the peak to peak particle displacement ampli¬ 
tude exceeds the length of the heat exchangers. Although to the 
knowledge of the author there is no clear experimental evidence 
for this transition effect in literature, some analogue findings are 
reported in the numerical studies of Worlikar and Knio [12] and 
of Marx and Blanc-Bennon [24]. The last authors, in particular, ob¬ 
served the same transition in behavior in the acoustically stimu¬ 
lated hydrodynamic energy flux flowing within the gas pores of 
an isolated stack (without heat exchangers). It is interesting to note 
as the effect (attributed by the authors to thermal wave anharmo- 
nicity) appeared at the drive ratio for which the peak-to-peak par¬ 
ticle displacement amplitude reached half the plate length. This 
condition could be retained as equivalent to the analogue put for¬ 
ward here for the HXs when considering that the plate edges of an 
isolated stack act as real HX fins and that at the operation point of 
the transition the stack may be considered as constituted by two 
HX fins (each of length L s /2) set in close contiguity. 

As reported above, the analogous approaching of a linear behav¬ 
ior is not as evident in q H as in q c . This fact could be explained 



Fig. 9. The heat loads per unit area on the cold and hot HX as a function of DR. In the 
insert their difference is plotted. Dotted and continuous lines result respectively 
from a quadratic (low DR) and linear (high DR) fit of the simulation points. 


observing that the hot HX have to sustain, in addition to the heat 
flux thermoacoustically pumped within the stack, the heat gener¬ 
ated by viscous dissipation over all the solid surfaces of the chan¬ 
nel (both stack and HXs). This contribute, which accounts for the 
larger heat transfer area of the hot HX compared to the cold one, 
has a well-known quadratic dependence on DR and should corre¬ 
spond in our simplified model to the difference q H - q c . This last 
is plotted in the insert of Fig. 9 where the quadratic dependence 
on DR is clearly observable. 

3.3. Effect ofHX length and fin interspacing 

The HX length along the axial direction is a critical parameter in 
HX design because it determines the available heat transfer area 
and also accounts for the rate of viscous dissipation. To study the 
effect of the cold HX length L c on the cooling load Q c , simulations 
are performed at fixed stack position and blockage ratio and at 
varying DR, y 0 , L c and L H values. A typical result is shown in 
Fig. 10 where Q c is plotted as a function of L c at selected DR values 
(runs 15-17). Simulations reveal that Qc increases quickly in the 
range 0 < L c < Xj becoming almost constant for L c >x A . At the con¬ 
sidered operating conditions and for L c l2xr ~ 1 the cooling load 
amounts at about 93% of the maximum cooling load evaluated 
by calculating the total longitudinal enthalpy flux at the midpoint 
of the stack (a choice justified by the fact that the stack plate is 
longer than the particle displacement for all the DR values consid¬ 
ered). As pointed out by Swift [1], this should be the optimal con¬ 
dition for efficient inter-element heat transfer. This conclusion is 
based on the argument that a particle can transport heat over a dis¬ 
tance not exceeding 2xr and that the total amount of particles able 
to participate to the inter-element heat transfer corresponds to the 
gas portion comprised within a distance 2x^ from the stack edge. 
We remarkably note as for all the considered drive ratios halving 
the plate length starting from the above invoked optimal value 
(2x, -> X!) produces a reduction of Q c of less than 1% (from about 
93% to 92%). This result indicates, as already highlighted by Hofler 
[25], that the fin length could be chosen substantially shorter than 
the peak-to-peak particle displacement amplitude without com¬ 
promising the HX performance and with great benefits in terms 
of viscous losses reduction. An experimental evidence for this find¬ 
ing can be found in [23] where it was recognized that doubling the 
length of the cold HX of a TA engine did not enhanced perfor¬ 
mance, even when the displacement amplitude much exceeded 
the length of the single HX. 



Fig. 10. The cooling load as a function of the length of the cold HX at selected drive 
ratios. Continuous lines are guides for the eye. 
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The dependence of the cooling load Q c on the plate/fin spacing 
2y 0 at selected L c values is shown in Fig. 11 (runs 18-23). At a fixed 
y 0 , the curve magnitude grows at increasing the fin length but the 
growth rate reduces with L c until the curves are being overlapping 
each other for high enough L c values. This trend is consistent with 
the results shown in Fig. 10 and should correspond to the flat por¬ 
tion of the curves attained for very long fins ( L c » Xi). Most impor¬ 
tant is that all curves peak at a well-defined value of the fin 
interspacing. The peak location results to be a slightly increasing 
(linear) function of the fin length, as evidenced in Fig. 11 by the 
straight line passing through the locus of maxima. This should cor¬ 
respond to a reduction of the effective length available for heat 
transfer (i.e. of the HX length) for short plate spacing (y 0 /<5 K < 1). 
The optimum fin/plate spacing ranges between 2y 0 = 2.98i5 K (for 
L c l 2x! = 0.079) and 2y 0 = 3.345 K (for Lc/2x-i = 1.9) for 12 = 0.76. 
These results agree with the general design procedure of selecting 
the plate spacing of the stack of thermoacoustic refrigerators in the 
interval 25 K ^ 2y 0 < 4 S K [26]. 

The difference between the cooling load and the total axial en¬ 
thalpy flux evaluated at the midpoint of the stack (maximum cool¬ 
ing load) reveals that the HX performance is affected by some kind 
of thermal losses. To get insight into this issue the transverse com¬ 
ponent of the heat flux density at the gas-solid interface (y = y 0 ) is 
plotted in Fig. 12 as a function of the axial coordinate (runs 24, 25). 
In both the upper (L c l2x^ = 0.53) and the lower (L c /2x! =2.25) 
curve illustrated the heat flux exhibits a sharp peak near the fin 
ends facing the duct (outer ends). These are therefore the main 
areas where a net heat exchange between the fluid and the solid 
takes place The lower plot evidences as the heat flux density in¬ 
verts its sign near the fin ends facing the stack (inner ends) of both 
cold and hot HX giving rise to a reverse flux which decreases the 
net heat transfer between the gas and the HXs, thus worsening 
the HX performance. Smaller transverse heat fluxes, qualitatively 
similar to those observed in thermoacoustic couples [27], are also 
found near the plate ends but they quickly vanish in the middle 
plate regions. The extent of the reverse heat fluxes, which are qual¬ 
itatively similar to those found by Besnoin and Knio [13,14], re¬ 
duces with DR till to vanish at high drive ratios when L c / 2xi -c 1, 
as evidenced by the upper curve of Fig. 12. An analogous trend is 
observed if DR is held fixed and L c is progressively reduced. The re¬ 
verse flux effect or, in general, the poor efficiency of the fin areas 
facing the stack (inner ends) in transferring heat could account 



Fig. 11. The cooling load as a function of the plate/fin interspacing at selected L c 
values and for DR = 4.93. Continuous lines result from the fit of the simulation 
points by an arbitrary peak function. The straight line passes through the locus of 
maxima. 



Fig. 12. The time-averaged transverse heat flux density at plate/fin surface (y = y 0 ) 
as a function of the axial coordinate for two DR values. 

for the previously discussed high performance of HXs whose 
length is considerably smaller than 2xj. Reducing the length of 
the HX, in fact, results in decreasing the fin areas of poor efficiency 
(inner ends) so, even if the available heat transfer area becomes 
smaller, the heat transferred per unit area is increased. This is 
clearly observable in Fig 13 where the mean heat flux per unit sur¬ 
face area | 'q^/lA exchanged by the cold fin with the gas is plotted 
as a function of the fin length at selected drive ratios (runs 15-17). 
The area over which the HX is affected by the reverse flux also var¬ 
ies with the fin/piate spacing as illustrated in Fig. 14 (runs 18, 19, 
21) where the longitudinal profile of the transverse heat flux den¬ 
sity at the gas-solid interface (y =y 0 ) of the cold HX fin is plotted at 
selected y 0 values. The comparative analysis of this figure with 
Fig. 11 informs that the points where Q c peaks are almost coinci¬ 
dent with the ones for which the reverse flux area results mini¬ 
mized. The lower plot of Fig. 12 reveals however that even when 
the reverse flux is eliminated, minor transverse heat fluxes local¬ 
ized at the stack edges persist. 

The thermal losses above discussed affecting the HX-stack junc¬ 
tion seems therefore to constitute a permanent feature of the inter¬ 
element heat transfer process and cannot be totally avoided. The 
results of the present study informs however that their strength 
can be minimized by suitable choices of the fin length and spacing 
which thus have to be simultaneously optimized. 



Fig. 13. The mean heat flux per unit surface area exchanged by the cold fin with the 
gas at selected drive ratios. 






















4528 


A. Piccolo / International Journal of Heat and Mass Transfer 54 (2011) 4518-4530 



Fig. 14. The longitudinal profile of the time averaged transverse heat flux density at 
the solid surface of the cold fin at selected fin spacing. 


3.4. Effect of blockage ratio 

The analysis performed in the last section indicates that the 
optimum value of plate/fin spacing for which the cooling load 
peaks is around 3<5 K . This result was derived considering a constant 
blockage ratio, that is, in each run the plate thickness was “ad¬ 
justed” to give rise for the selected y 0 to the constant value 
Q = 0.76. The analysis, therefore, did not take into account the ef¬ 
fect of varying plate/fin spacing on particle axial velocity and effec¬ 
tive cross sectional area available for heat transport. In order to 
estimate the effect of blocked ratio on the total cooling load simu¬ 
lations are performed at varying y 0 for a fixed value of the plate 
thickness / (roughly equal to the thermal penetration depth in 
the fin material: 0.0005 m). Also, the total number of channels 
comprised in a cross section of the stack are calculated assuming 
that this last fill in all the resonator cross-section, taken for sim¬ 
plicity quadratic in shape and of unitary surface area. The total 
cooling load sustained by the stack is then obtained by multiplying 
2Q c for the resulting number of channels. Results are shown in 
Fig. 15 where the total cooling load is reported as a function of 
y 0 (and Q) at selected L c values (runs 26-28). All curves result 
much more peaked compared to the ones reported in Fig. 11 and 
the peak location (also in this case slightly dependent on L c ) falls 
at lower y 0 values, i.e., at a fin/plate interspacing roughly equal 
to 2S k (to which the optimal value Q 0.59 corresponds). The dif¬ 
ferent behavior can be interpreted observing that the lower the 
plate interspacing the higher the number of channels of the stack. 
This compensates for the decreasing of the heat flux pumped in a 
single channel which is observed in Fig. 11 to take place for 
y 0 < 1.5<5 k . Simulations reveal that a change in plate thickness mod¬ 
ifies the peak height and associated optimum blockage ratio but 
not the peak location. For example, doubling the plate thickness 
(for Lc/2xi = 1) reduces the peak height of about 9% and the optimal 
Q by about 29% leaving the peak location unchanged (results not 
shown). Based on the results of this section the optimal plate/fin 
interspacing in stack-based thermoacoustic device should fall in 
the narrower range 2S K < 2y 0 ^ 3<5 K . 

3.5. Derivation of heat transfer correlations 

In this section simulations are performed with the intent to de¬ 
velop heat transfer correlations which could be applied for the esti¬ 
mation of the gas-side heat transfer coefficient h at the HX walls. At 
the same time, the comparison of the theoretical predictions 
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Fig. 15. The total cooling load per unit cross sectional area of the stack as a function 
of blockage ratio at selected lengths of the cold HX. 


against experimental data found in literature is considered as 
model validation. 

Due to the existence of reverse heat flux localized at the inner 
edges of the HX fins the area-averaged h values over the cold HX 
are computed as indicated in [13] by the formula 


h = 


(Qy) 

(Thxc) - (T g c) 


(31) 


where 


1 f Lc 

(<? y) = r / 9y(X,y o)dx 
Lc Jo 


(32) 


Analogous formulas have been applied to the hot HX. The h data 
thus obtained are converted to non-dimensional Nusselt numbers 
based on the hydraulic diameter of the pore D h (=4y 0 ) and then to 
Colburn-j factors 


NUn = 


hD„ 


Jc,D l 


Nu n 


R e,. D Pr 1/3 


(33) 


Re] D [=t'xi(0)D h /v] being the acoustic Reynolds number defined in 
terms of D h and of the amplitude of the acoustic velocity evaluated 
at gas centerline. The results of about 40 simulations (runs 29-37) 
carried out at varying DR, L c , L H , y 0 , and x s values are illustrated in 
Fig. 16 where the resulting Colburn-j factors are plotted as a func¬ 
tion of Re] D . In the same graph the experimental measurements 
of Nsofor et al. [28], Peak et al. [29], Mozurkewich [9] and Brewster 
et al. [22] are reported. All HXs involved in these works are of the 
finned-tube type (except the one in [9] consisting of a copper screen 
soldered to tubes referred to by the author as “configuration C”) so 
the measured h values therein reported are considered suitable to 
make comparison with the predictions of the proposed model. Pre¬ 
liminary, in order to make the different set of data consistent, the 
experimental Nusselt and Reynolds numbers have been re-scaled 
by the same length parameter (4y 0 ). In addition, the rms velocity 
amplitudes entering in the Reynolds numbers reported in [5] have 
been “corrected" by the factor V2 and by the blockage ratio £2. 

The mean deviation of the predicted j c values from the experi¬ 
mental ones amounts to 49% for the data of Peak et al„ 36% for the 
data o Brewster et al., 36% for the data of Mozurkewich and 56% 
for the data of Nsofor et al. Deviations in the moderate and high Rey¬ 
nolds number regime (Re] > 500) are expected because of the in¬ 
creased weight of entrance/exit effects at the resonator-HX cross 
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section interface. These effects are responsible of complex non-lin¬ 
ear temperature and flow patterns (turbulent and vorticity flows 
[6].) which highly affect heat transfer rates causing, among other 
things, minor losses and heat leakage to the duct. The high deviation 
found at low Rei numbers (Re! < 500) between simulations and 
experimental results of Nsofor et al. is quite surprising because 
the linear theory, on which the present model relies upon, is gener¬ 
ally observed to match experiments quite well at low velocity 
amplitudes. In the evaluation of these results it has to be pointed 
out however that, although proper scaling procedures have been ap¬ 
plied, the operating and geometrical conditions specified in the sim¬ 
ulations are different from those of the selected experiments. It’s 
remarkable to note, anyway, that al low Rei values (Re! < 200) the 
simulation points match exactly the correlation curve derived by 
Peak et al. from the fit of their data. On the same curve, on the other 
hand, the data of Mozurkewich tend to collapse at decreasing Rei. 

4. Conclusions 

The performance of parallel-fin thermoacoustic HXs coupled to 
a stack in a simplified configuration is investigated through a 2D 
numerical model based on the classical linear thermoacoustic the¬ 
ory. Details of the transverse heat transfer at the HXs are high¬ 
lighted and optimization information for their design are derived. 
The following conclusion can be drawn: 

1. The magnitude of the heat transfer coefficient from the second¬ 
ary fluid-side (If) should be fully taken into account because it 
affects cooling load through its influence on the temperature 
distribution across the stack. High U values are desirable to 
improve the performance of the device. 

2. Fin length along the axial direction of particle oscillation can be 
chosen considerably lower than the peak-to-peak displacement 
amplitude without compromising the HX performance and 
with great benefits in terms of viscous losses reduction. 

3. Thermal losses (reverse heat fluxes) localized at the stack-HX 
junctions degrade the HXs performance reducing the useful 
cooling load. These effects could account for about 7% of the 
deviations found between predictions of the linear theory and 
experimental measurements. Suitable choices of the fin spacing 
and fin length can contribute to minimize this detriment effects 
enhancing the effectiveness of the HXs. 

4. Optimum fin/plate interspacing should fall in the range 
23 k ^ 2y 0 ^ 4S k . Selection of the best value is case dependent 
and should also take into account the simultaneous effect of 
blockage ratio, fin length and fin thickness. 


5. Heat transfer coefficients from the gas-side can be predicted 
with a confidence of 36% and 49% respectively at moderate 
and high acoustic Reynolds numbers. The high deviation found 
at low Rei with available experimental data requires further 
investigation. 

The simplified computational method proposed in this work 
suffers by several limitations. In its current form the model is re¬ 
stricted to the acoustic and short stack approximations. Further¬ 
more, the HXs are modeled as coupled to the stack in a 
simplified configuration while their arrangement in real TA devices 
can be more complex. The heat transport equations implemented 
in the code relies on the classical form of the function Ti whose 
accuracy is questionable near the pore terminations. The model 
does not take into account for entrance and exit effects which 
may become relevant near sudden changes of cross section like 
those existing at the ends of the HXs. Further analysis is needed, 
particularly of the fluid and energy flows near the HX pore ends, 
for overcoming these shortcomings. 
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